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1.  INTRODUCTION 


Are  the  reactions  that  occur  behind  the  shock  front  in  a  detonation  due  to  rapid  heating  of  molecules 
to  dissociation  after  shock  wave  passage,  or  are  the  bonds  mechanically  broken  after  shock  compression? 
Or  do  reactions  result  from  a  combination  of  both  processes?  What  are  the  consequences  for  the 
properties  of  the  detonation  wave  of  each  of  these  mechanisms?  The  mechanism  of  detonation  remains 
unknown,  despite  numerous  investigations  over  the  past  century  (Pickett  and  Davis  1979;  Pickett  1985). 
The  major  obstacle  in  determining  reaction  mechanisms  from  experiment  lies  in  the  short  time  and  spatial 
scales  over  which  detonation  occurs  and  is  complicated  by  the  extreme  pressure  and  energy  releases 
accompanying  the  detonation. 

Recent  advances  in  spectroscopic  methods  are  beginning  to  allow  the  direct  measurements  of 
molecular  changes  in  a  system  that  is  shocked,  toward  the  goal  of  determining  the  atomistic  processes 
involved  in  the  detonation  (Gupta  1995).  Others  are  investigating  energy  transfer  rates  (Pried  and 
Ruggiero  1994)  and  mechanisms  (Dlott  and  Payer  1990;  Tokmakoff,  Payer,  and  Dlott  1993;  Chen,  Tolbert, 
and  Dlott  1994;  Chen,  Hong,  Hill,  and  Dlott  1995)  in  energetic  materials  and  possible  relationships  to  the 
initiation  of  the  detonation.  These  studies  are  providing  critical  pieces  to  a  multifaceted  chemical  problem 
that  is  not  completely  understood.  The  measurements  of  molecular  changes  of  explosives,  however,  have 
not  been  made  under  the  extreme  conditions  associated  with  detonation.  Direct  measurement  of  chemical 
reactions  occurring  during  this  destructive  and  rapid  event  remains  a  formidable  hurdle  that  has  not  been 
overcome. 

As  the  development  of  diagnostic  tools  to  measure  these  ultrafast  events  has  progressed,  computational 
resources  and  theoretical  methods  to  model  atomistic  processes  involved  in  tfie  detonation  have  advanced 
as  well.  Early  reactive  models  used  in  molecular  dynamics  simulations  described  molecular  crystals 
consisting  of  diatomics  in  metastable  bound  states  that,  when  shocked,  dissociated  exothermically  to  form 
the  mote  stable  monatomic  products  (Karo,  Hardy,  and  Walker  1978;  Tsai  and  Trevino  1981).  A  more 
recent  study  (Kawakatsu  and  Ueda  1988, 1989;  Kawakatsu,  Matsuda,  and  Ueda  1988),  in  which  molecular 
dynamics  simulations  of  the  detonation  were  compared  to  hydrodynamic  predictions,  used  an  isomerization 
of  the  reactant  as  the  energy-release  reaction  that  drives  the  detonation.  Although  these  models  predicted 
reasonable  pressure  and  temperature  rises  when  compared  to  observed  data  from  actual  detonations,  the 
chemistry  that  drives  the  detonation  in  these  models  is  qualitatively  unrealistic.  A  next  step  toward  a  mote 
realistic  description  of  the  chemical  processes  occurring  in  a  detonation  would  be  one  in  which  the  energy 
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release  is  due  to  endothermic  dissociation  of  the  reactants  followed  by  exothermic  association  of  the 
fragments  to  form  more  stable  products.  One  must  proceed  in  this  process  slowly,  that  is,  introducing  a 
minimum  of  additional  complexity  in  the  reaction  mechanisms  so  as  not  to  mask  incremental  information 
which  might  result. 

Brermer  and  coworkers  (Brenner  1992;  Robertson  et  al.  1992;  White  et  al.  1992, 1994;  Brermer  et  al. 

1993) ,  in  attempts  to  describe  mote  realistically  exothermic  chemical  reactions  occurring  in  a  shock  wave, 
developed  a  series  of  models  based  on  gas-phase  reactivity  schemes.  The  models  describe  a  heteronuclear 
diatomic  molecular  crystal  that  requires  energy  to  break  the  molecular  bonds  in  the  low-pressure  crystal, 
but  releases  substantial  energy  upon  formation  of  homonuclear  products.  Two-dimensional  molecular 
dynamics  simulations  using  their  first  model  (Brenner  1992;  Robertson  et  al.  1992;  White  et  al.  1992, 

1994)  (denoted  Model  0  hereafter)  appeared  to  sustain  shock  waves  that  were  driven  by  exothermic 
reactions.  However,  the  macroscopic  property  profiles  (pressure,  density,  and  temperature)  of  the 
detonating  system  were  flat-topped  and  split  (Brermer  1992;  Robertson  et  al.  1992;  White  et  al.  1992, 
1994) — ^this  behavior  has  not  been  observed  in  experiment  The  results  of  a  second  study  (Brermer  et  al. 
1993)  using  a  similar  model  (the  only  difference  being  in  the  molecular-size  parameter)  showed 
macroscopic  shock  profiles  that  agreed  well  with  the  simplest  model  of  detonation  (Pickett  and  Davis 
1979;  Pickett  1985),  in  which  reaction  takes  place  immediately  behind  the  shock  fiont  We  will  denote 
this  model  hereafter  as  Model  I.  Peatures  of  these  models  win  be  discussed  below. 

The  studies  described  in  this  paper  and  an  accompanying  paper  (Rice  et  al.,  in  press)  were  partially 
prompted  after  our  attempts  to  reproduce  simulation  results  (published  in  Brermer  et  al.  1993)  using  Model 
I  failed.  Additionally,  Model  I  (as  well  as  Model  0)  had  features  in  the  potential  energy  surface  that 
seemed  undesirable  to  us,  as  described  later  in  this  report.  We  made  changes  to  Model  I  by  removing 
these  undesired  features  (the  revised  model  denoted  as  Model  n  hereafter)  and  then  performed 
two-dimensional  molecular  dynamics  simulations  of  the  detonation  of  this  model  energetic  crystal  initially 
at  low  temperature,  low  pressure.  The  equation  of  state  and  Hugoniot  curve  for  this  model  were  also 
obtained  from  two-dimensional  molecular  dynamics  simulations  at  appropriate  temperatures  and  pressures. 
These  were  used  to  make  hydrodynamic  predictions  of  detonation  pressures,  densities,  and  velocities.  The 
hydrodynamic  predictions  were  in  good  agreement  with  the  results  of  our  simulations  of  the  detonation, 
indicating  that  our  model  and  the  method  of  molecular  dynamics  are  consistent  with  the  conservation  laws 
of  the  process  of  detonation. 
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The  focus  of  the  study  reported  here  is  to  examine  the  microscopic  details  of  the  detonation  process 
and  extract  infonnation  concerning  the  mechanism  of  the  detonation.  We  will  first  provide  a  detailed 
analysis  of  the  interaction  potentials  for  Models  0, 1,  and  II  and  show  that  the  reaction  mechanisms  of 
these  systems,  when  subjected  to  shock,  are  strongly  dependent  on  the  functional  form  of  the  models. 
Then,  using  the  results  from  two-dimensional  molecular  dynamics  simulations  reported  in  our 
accompanying  paper  (Rice  et  al.,  in  press),  we  will  dissect  the  reaction  zone  and  follow  the  changes  in 
molecular  properties  throughout  this  region.  This  will  reveal  the  mechanism  of  detonation  for  this  model. 

2.  MODELS 

2.1  Models  0  and  I.  Models  0  and  1,  used  by  Brenner  (Brenner  1992;  Robertson  et  al.  1992;  White 
et  al.  1992;  White  et  al.  1994;  Brenner  et  al.  1993),  are  two-dimensional  crystals  of  diatomic  molecules 
arranged  in  a  herringbone  lattice.  The  interaction  potential  that  is  used  to  describe  these  crystals  is 

^  k  ('ii)  [''r  (■'ij)  -  ''a  (■■ij)  ]  *  vJP*  (ry)  j  ,  (1) 

1  J>1 


the  sums  being  over  the  N  atoms  comprising  the  model.  The  functional  forms  and  parameters  for  the 
terms  in  Equation  (1)  are  in  given  in  Table  1.  The  expression  in  Equation  (1)  is  based  on  ideas  proposed 
and  developed  by  AbeU  (1985)  and  Tersoff  (1988).  The  left-most  terms  within  the  suimnations  in 
Equation  (1)  make  up  the  intramolecular  interaction  potential  of  diatomic  molecules,  and  the  second  term 
within  the  sums,  is  an  intermolecular  interaction  term.  The  function  f^(r^)  smoothly  attenuates  the 
molecular-bonding  interaction  to  zero  at  3.0  A.  The  bond-order  function  B^,  which  ranges  in  value  from 
zero  to  one  depending  on  the  local  atomic  enviromnent,  introduces  many  body  effects  by  modifying  the 
attractive  term  of  the  molecular-bonding  portion  of  Equation  (1),  Vy^(rjj).  For  an  isolated  diatomic,  the 
value  of  is  one,  and  Vy^(ry)  will  have  a  maximum  contribution  to  Equation  (1)  for  the  ij  pmr  at  the 

specified  intemuclear  distance.  If  the  diatomic  is  closely  surrounded  by  many  other  atoms,  as  in  the  high- 
dense  region  of  a  shock-compressed  crystal,  the  value  of  B^  decreases,  depending  on  the  number  and 
location  of  the  nearest  neighbors  surrounding  the  ij  pair.  This  decrease  in  B^  will  correspondingly 
decrease  the  attractive  interaction  for  the  ij  pair.  This  mbric  is  intended  to  mimic  modifications  of  the 
bonding  character  upon  increased  density.  As  is  made  clear  by  the  notation  for  the  parameters  given  in 
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Table  1.  Parameters  and  Functional  Forms  Used  for  the  Potential  Energy  Expressions  in  Equations  (1)  and  (2) 
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Reference  Brenner  (1992).  There  is  a  typographical  error  in  Table  1  of  this  reference  for  the  By  definition  (Brenner,  private  communication).  It  is  corrected  in  Brenner 
et  al.  (1993). 

Brenner  et  al.  (1993). 


Table  1,  this  fonn  of  the  potential  is  capable  of  describing  diatomic  molecules  fonned  by  atoms  of  types 
A-B,  A- A,  and  B-B  as  well  as  intermolecular  interactions  between  these  various  atom  types.  The  curves 
in  Rgure  1(a)  illustrate  the  effect  of  the  term  on  the  intramolecular  term  for  an  A-B  pair.  The  1 1 
curves  in  this  figure  correspond  to  the  A-B  intramolecular  interaction  with  values  ranging  from  0  to  1. 
This  figure  shows  that  the  ij  attraction  is  decreased  by  over  one-half  for  B^  =  0.7,  and  the  molecular 
interaction  appears  to  be  completely  repulsive  for  B^<0.1.  Also,  the  positions  of  the  potential  energy 
minima  increase  with  decreasing  B^.  The  behavior  of  this  function  with  decreasing  B^  indicates  that  as 
the  local  atomic  environment  becomes  increasingly  dense  aroimd  an  ij  pair,  the  i  and  j  atoms  will  be 
"forced  apart"  by  the  increased  repulsion  that  can  be  attributed  to  this  functional  form. 

Figure  1(b)  shows  the  effect  of  the  B^  term  on  the  intramolecular  term  for  an  A-A  (or  B-B)  pair. 

(The  masses  and  homonuclear  interactions  of  A  and  B  atoms  are  the  same  in  Model  0).  As  in  Figure  1(a), 
the  intramolecular  attraction  is  decreased  by  over  one-half  for  B^ = 0.7  and  becomes  completely  repulsive 

for  B^<0.1.  However,  due  to  the  differences  in  weU  depths  for  the  intramolecular  potential  for  A-A  (or 
B-B)  pairs  and  the  A-B  pairs,  the  attractive  interactions  for  A-A  pairs  are  considerably  greater  than  A-B 
pairs  for  most  values  of  B^. 

Brenner  and  cowoikers  found  that  Model  0  produced  "flat-topped,  split  shock-wave  structure"  of  the 
detonating  crystal  that  is  due  to  a  "high-pressure,  dissociative  phase  transition  unintentionally  introduced 
through  our  initial  parameterizations"  (Brenner  1993).  In  our  preliminary  studies  using  Model  0,  we  were 
able  to  reproduce  each  feature  of  the  potential  energy  surface  that  was  reported  (Brenner  1992;  Robertson 
et  al.  1992;  White  et  al.  1992;  White  et  al.  1994),  including  sound  speed,  equilibrium  nearest-neighbor 
distance,  barrier  to  collinear  reaction  of  A  +  A-B,  and  shock  wave  velocities.  We  also  reproduced  the  flat- 
topped  split  shock- wave  structures  alluded  to  in  footnote  13  of  Brermer  (1993). 

Brenner  et  al.  subsequently  published  results  fi'om  simulations  using  Model  I  in  which  the  "high- 
pressure,  dissociative  phase  transition"  does  not  occur,  and  the  detonation  profiles  show  a  single  shock 
front  that  is  followed  by  a  reacting  flow  and  rarefaction  wave  (Brenner  1993).  The  results  using  Model  I 
were  those  that  we  could  not  reproduce. 

Model  I,  as  repotted  in  Brenner  (1993),  is  the  same  as  Model  0  with  one  difference;  that  is,  the  value 
of  tg.  In  Model  0,  r^  had  a  value  of  1.0  A;  r^  in  Model  I  is  1.2  A.  Thus,  the  changes  in  intramolecular 
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Figure  1.  Intramolecular  interaction  potentials  of  Model  0  for  different  values  (Equation  1)  as 
functions  of  intemnclear  distance  for  (a)  A-B  interactions:  (b)  A-A  interactions.  Magnitude 
of  the  values  are  denoted  in  the  legend. 

6 


HUH 


interaction  potential  with  for  Model  I  are  similar  to  those  of  Model  0,  shown  in  Figure  1,  with  the 
equilibrium  distances  shifted  by  0.2  A  for  each  curve. 

Using  the  parameters  of  Model  I,  we  have  performed  molecular  dynamics  simulations  in  an  attempt 
to  reproduce  results  of  Brenner  (1993).  However,  we  once  again  saw  the  flat-topped,  split  shock-wave 
structure  that  was  seen  in  Model  0.  Figure  2  is  a  snapshot  of  the  detonating  system  at  15  ps  resulting 
ftom  our  simulations  using  Model  I.  This  snapshot  shows  that  four  distinct  regions  exist:  the  undisturbed 
crystal,  a  compressed  crystal  in  which  the  molecules  are  reoriented,  but  do  not  react;  the  "high-pressure, 
dissociative  phase  transition"  (Brenner  1993)  in  which  molecular  identity  is  lost,  and  the  rarefaction  region, 
consisting  of  vibrationaUy  excited  homonuclear  products.  The  velocities  of  the  two  leading  discontinuities 
are  8.1  and  6.8  km/s,  respectively,  whereas  Brenner  (1993)  reported  the  velocity  of  the  single  shock  fiont 
as  9.3  km/s. 

We  do  not  know  the  source  of  the  discrepancies  between  our  calculations  and  those  of  Brermer  (1993). 
Several  details  about  the  calculations,  particularly  regarding  the  initial  conditions,  were  not  given  in 
Brenner  (1993),  including  the  equilibrium  molecular  orientation,  lattice  constants,  temperature  of  flie 
undisturbed  crystal,  and  number  of  atoms  in  the  system.  We  therefore  had  to  attempt  to  reproduce  tiie 
simulations  reported  in  Brenner  (1993)  with  an  incomplete  set  of  information.  We  determined  the 
equilibrium  low-pressure  (1  atm),  low  temperature  (20  K)  crystal  structure  for  Models  0  and  I  using  NPT 
Metropolis  Monte  Carlo  sampling  as  described  in  Rice  et  al.  (in  press).  The  low-temperature,  ambient- 
pressure  lattice  parameters  of  the  crystal  defined  by  Model  I  in  the  x-  and  y  directions  are  a  =  4.16  A  and 
b  =  6.25  A,  respectively.  The  center-of-bond  (COB)  firactionals  for  the  two  molecules  (1  and  2)  in  the 
unit  cell  are  at  (0.25, 0.25)  and  (0.75, 0.75),  respectively.  The  equilibrium  A-B  bond  length  is  1.19945  A, 
and  the  molecular  orientation  of  molecules  1  and  2  are  30.51°  and  -30.51°,  respectively,  relative  to  the 
crystal  x-axis.  The  details  of  the  molecular  dynamics  simulations  are  given  in  Section  3.4  of  Rice  et  al. 
(in  press). 

Upon  finding  these  discrepancies  in  our  calculations  when  compared  to  those  published  in  Brenner 
(1993),  we  attempted  to  understand  the  source  of  the  "dissociative  jAiase  transition,"  which  seemed  to 
cause  the  problem  of  the  flat-topped  shock  stracture.  We  began  our  analyses  by  examining  the  three-body 
molecular  interaction  term  that  includes  for  Model  I.  We  generated  contour  plots,  shown  in  Rgure  3 

of  the  collinear  reactions: 
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contours  (in  eV)  for  the  following  collinear  three-bodv  reactions  usin 


A  +  BA  -->  A-B  +  A:  fc)  A  +  A-A  -->  A-A  +  A. 


Units  for  the  intemuclear  distances  along  both  axes  are  in  A. 
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A  +  A-B  ->  A-A  +  B 

AE  =  -3.0  eV 

(D 

A  +  B-A  ->  A-B  +  A 

AE  =  0.0  eV 

(n) 

A  +  A-A  ->  A-A  +  A 

AE  =  0.0  eV 

(HI) 

Note  that  because  the  masses  and  bond  energies  of  A-A  and  B-B  are  the  same  for  Model 
(I)-(ni)  have  the  same  potential  energy  features  as  the  following  collinear  reactions: 

I,  that  reactions 

B  +  B-A 

AE  =  -3.0  eV 

(IV) 

B  +  A-B 

AE  =  0.0  eV 

(V) 

B  +  B-B. 

AE  =  0.0  eV 

(VI) 

As  previously  stated,  we  were  able  to  reproduce  the  contour  plot  and  barrier  to  Reaction  (1)  for 
Model  0,  given  in  Brenner  (1992),  Robertson  et  al.  (1992),  and  White  et  al.  (1992,  1994).  Figure  3(a), 
the  contour  of  Reaction  (I)  for  Model  I,  is  similar  in  feature  to  the  contour  of  Model  0  given  in  Brenner 
(1992),  Robertson  et  al.  (1992),  and  White  et  al.  (1992,  1994);  we  calculate  a  barrier  to  reaction  upon 
collinear  j^proach  as  0.10  eV,  whereas  the  value  reported  in  Brenner  (1993)  is  0.08  eV  (this  may  not  be 
significant). 

However,  as  is  evident  in  Figures  3(b)  and  3(c),  the  trimers  ABA  and  AAA  exist  upon  linear 
approach.  (Ehie  to  the  features  of  the  potential,  this  indicates  that  the  trimers  BAB  and  BBB  also  exist) 
We  characterized  these  species  through  normal  mode  analyses,  after  finding  the  equilibrium  structures  of 
these  trimers  using  the  Newton-Raphson  energy  minimization  method  (Margenau  and  Murphy  1943).  The 
geometries  and  vibrational  frequencies  of  these  species  are  given  in  Table  2.  Also,  we  characterized  the 
structures  corresponding  to  the  saddle  points  apparent  in  these  figures.  We  were  able  to  reproduce  the 
features  of  the  A-A  and  A-B  diatomics  that  were  reported  in  Brenner  (1993),  as  well  as  the  crystalline 
binding  energy  of  0.04  eV  per  molecule  and  nearest-neighbor  distances  of  3.3  A.  We  were  unable  to 
locate  any  nonlinear  trimers  AAA  or  ABA  or  any  of  the  trimers  AAB  or  BBA  using  Equation  (1)  for 
either  Models  0  or  1.  Higher  order  multimers  were  not  found. 

The  relative  energy  features  of  the  potential  energy  surface  for  Model  I  are  illustrated  in  Figure  4. 
These  schematics  show  relative  energies  of  the  trimers  and  the  saddle  points  leading  to  their 
formation/decomposition.  The  shallow  minima  of  the  trimers  relative  to  the  saddle  points  leading  to 
dissociation  suggest  that  the  trimers  would  be  short-lived  in  a  high-energy  environment  such  as  that  which 
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Table  2.  Geometric  Parameters  and  Harmonic  Vibrational  Frequencies  of  Critical  Points 
on  the  Potential  Energy  Surfaces  for  Models  0, 1,  and  n 


Species 

Intemuclear  Distances 

(A) 

Harmonic  Frequencies 
(cm'b 

Energy 

(eV) 

Model  0 

A-B 

R(A-B)  =  l.(X)000 

1064 

-2.0 

A-A  (B-B) 

R(A-A)  =  1.00000 

1682 

-5.0 

A-A-A  (B-B-B) 

R(A-A)  =  1.24082 

200,  296,  739 

-3.94069 

A-B-A  (B-A-B) 

R(A-B)  =  1.23822 

150,  184,  466 

-1.58114 

A  +  A-B  “>A-A  +  B 

R(A-A)  =  2.15000 
R(A-B)  =  1.01298 

217j,  982 

-1.89457 

A  +  B-A  ~>A-B-A 

R(A-B)  =  1.37625 
R(B-A)  =  1.14901 

158i,  171,  570 

-1.57884 

A  +  A-A  ->A-A-A 

R(A-A)  =  1.35566 
R(A-A)  =  1.16141 

225i,  280,  866 

-3.93764 

Model  I 

A-B 

R(A-B)  =  1.2 

1064 

-2.0 

A-A  (B-B) 

R(A-A)  =  1.2 

1683 

-5.0 

A-A-A  (B-B-B) 

R(A-A)  =  1.38609 

195,  690,  863 

-4.42361 

A-B-A  (B-A-B) 

R(A-B)  =  1.38554 

122,  439,  545 

-1.77019 

A  +  A-B  ->A-A  +  B 

R(A-A)  =  2.28375 
R(A-B)  =  1.21310 

257r,  974 

-1.90225 

A  +  B-A  ~>A-B-A 

R(A-B)  =  1.87013 
R(B-A)  =  1.24572 

230i,  847 

-1.64694 

A  +  A-A  ->A-A-A 

R(A-A)  =  1.87013 
R(A-A)  =  1.24572 

365i,  1339 

-4.11734 

Model  n 

A-B 

R(A-B)  =  1.35 

727 

-1.0 

A-A 

R(A-A)  =  1.2 

1626 

-5.0 

B-B 

R(B-B)  =  1.5 

587 

-2.0 

A  +  B-A  “>A-B  +  A 

R(A-B)  =  1.64967 

300j,  275 

-0.57237 

B  +  A-B  ~>B-A  +  B 

R(A-B)  =  1.64967 

171j,  157 

-0.57237 

A-B  +  B  ->A  +  B-B 

R(A-B)  =  1.68460 
R(B-B)  =  1.79843 

367i,  365 

-0.82174 

B  +  B-B  ->B-B  +  B 

R(B-B)  =  1.79967 

242i,  222 

-1.14474 

A  +  A-A  — >A-A  +  A 

R(A-A)  =  1.49969 

670/,  631 

-2.86186 
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exists  behind  the  shock  front.  Additionally,  the  trimer  minima  will  be  affected  by  the  density,  since  trimer 
existence  is  due  to  the  intramolecular  interaction  term  in  Equation  (1).  The  effect  on  the  trimer  minima 
due  to  the  changes  in  density  as  the  shock  front  travels  through  the  region  wiU  affect  the  lifetime  of  the 
trimer  intermediate.  The  ability  of  Equation  (1)  to  form  trimers  leads  us  to  suspect  this  might  explain  the 
region  in  the  detonating  system  that  Brenner  et  al.  call  the  "dissociative  phase-transition"  (Biermer  1993). 
We  think  it  might  be  more  appropriate  to  label  that  region  as  an  associative  phrase-transition,  where  trimer 
formation  occurs,  leading  ro  subsequent  decomposition  to  the  more  stable  homonuclear  diatomic  species. 
In  any  case,  changing  the  value  of  the  parameter  r^  from  1.0  to  1.2  A  did  not  remove  the  source  of  the 
"flat-topped,  split  shock-wave  structure"  (Brenner  1993),  which  we  think  is  due  to  this  ability  to  form 
trimers. 

2.2  Model  n.  The  ability  of  trimer  formation  using  Models  0  and  I  indicated  to  us  that  saturation 
of  the  molecular  bond  was  not  correctly  accounted  for  in  Equation  (1)  with  those  sets  of  parameters.  It 
seems  appropriate  to  us  that  not  only  should  the  attractive  portion  of  the  potential  be  modified  according 
to  the  local  atomic  environment,  but  that  the  repulsive  wall  should  be  affected  as  well,  to  take  into  account 
saturation  of  tire  bonds.  We  have  therefore  modified  Equation  (1)  to 

V  -  E  E  [(2  -B;:)  Vj(r,j)  -  v*  (r,j)]  *  vj?*} ,  ® 


where  the  term  has  the  same  description  as  in  Equation  (1).  The  parameters  and  functional  form  for 
the  terms  in  Equation  (2)  are  given  in  Table  1.  We  foimd  two  discontinuities  in  the  weak  long-range  van 
der  Waals  interaction  for  Models  0  and  I;  one  at  r  =  2.91  A  and  the  other  at  r  =  7.32  A.  We  have 
corrected  the  discontinuity  at  r  =  2.91  A  by  determining  the  cubic  spline  coefficients  using  the 
requirements  that  at  r  =  1.75  A,  die  energy  and  energy  first  derivatives  must  equal  0.0,  and  at  r  =  2.91  A, 
the  energy  and  energy  first  derivatives  must  equal  that  produced  by  the  Lennard-Jones  potential  function. 
We  also  added  a  quintic  spline  to  the  Lennard-Jones  potential  in  the  region  of  7.31-7.32  A,  to  allow  a 
continuous  cutoff  of  die  intermolecular  interaction.  The  quintic  spline  coefficients  were  determined  in  the 
same  marmer  as  the  cubic  spline  coefficients  described  earlier,  however,  energy  second  derivatives  were 
required  to  match  at  the  functional  interfaces  at  7.31  A  and  7.32  A.  We  have  only  provided  four 
significant  figures  for  the  spline  coefficients  used  in  Equation  (2);  we  suggest  that  machine  accuracy  of 
these  coefficients  be  determined  by  the  interested  modeler. 


13 


We  have  also  changed  the  model  parameters  to  allow  different  exothermic  reactions,  and  set  the 
masses  of  particles  A  and  B  to  15  and  46  amu,  respectively,  in  order  to  introduce  mass  effects.  Also,  the 
molecular  sizes  of  the  A-B,  A-A,  and  B-B  diatomics  were  set  to  1.35, 1.2,  and  1.5  A,  respectively.  The 
reactions  that  this  model  describes  are: 


A-B  +  A  — >  A  +  B-A  AE  =  0.0  eV  (VB) 

B-A  +  A  — >  B  +  A-A  AE  =  -4.0  eV  (Vm) 

A-B  +  B  — >  A  +  B-B  AE  =  -1.0  eV  (DC) 

B-A  +  B  — >  B  +  A-B  AE  =  0.0  eV  (X) 

A-A  +  A  — >  A  +  A-A  AE  =  0.0  eV  (XI) 

B-B  +  B  — >  B  +  B-B  AE  =  0.0  eV  (XU) 


Geometric  parameters  and  vibrational  frequencies  of  the  stable  species  for  these  reactions  are  given  in 
Table  2.  We  also  show  contour  plots  of  the  coUinear  reactions  using  Equation  (2),  assuming  only  three 
atoms,  in  Figure  5.  We  did  not  find  trimer  or  other  multimer  fomiations  using  Equation  (2).  For  the 
coUinear  approach  of  A  to  A-B,  there  is  no  barrier  to  tq)proach.  The  transition  state  structures  and 
vibrational  frequencies  for  the  coUinear  reactions  are  also  given  in  Table  2. 

We  have  given  considerable  attention  to  describing  the  features  of  the  potential  energy  surface  for 
Model  n  assuming  oidy  three  atoms,  primarily  because  such  attention  was  given  in  die  previous  studies 
for  Models  0  and  I  (Brenner  1992;  Robertson  et  al.  1992;  White  et  al.  1992, 1994;  Brenner  1993)  and  we 
wished  to  provide  direct  comparison.  However,  in  the  high-dense  region  behind  the  shock  wave,  the 
potential  energy  interactions  (specificaUy  the  intramolecular  interactions)  wiU  be  affected  by  many  more 
atoms  than  three.  Therefore,  as  we  did  for  Model  0,  we  show  the  effect  of  the  intramolecular  interaction 
with  different  values  of  in  Figure  6,  which  correspond  to  varying  degrees  of  density  of  atoms 

surrounding  an  ij  pair. 

These  curves  differ  from  those  of  Model  0  (Figure  1)  in  that  the  weU  depths  are  decreased  by  one-half 
for  =  0.8  rather  than  B^  =  0.7.  AdditionaUy,  the  repulsive  interactions  for  these  B^  values  are  much 

greater  than  those  for  Model  0,  and  the  positions  of  the  minima  are  at  larger  distances  with  decreasingB^ 
than  those  for  Model  0.  This  model  is  analogous  to  Models  0  and  I  in  that  the  atoms  are  forced  apart  with 
decreasing  B^. 
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Figure  5.  Potential  energy  contours  (in  eV)  for  the  following  collinear  three-body  reactions  using 
Model  n:  (a)  A  +  A-B  ->  A-A  +  6:0)1  A  +  B-A  ->  A-B  +  A:  (cl  A  +  B-B  ->  A-B  +  B: 
(d)  B  +  A-B  — >  B-A  +  B:  (e)  A  +  A-A  — >  A-A  +  A;  (f)  B  +  B-B  — >  B-B  +  B.  Units  for 
the  intemuclear  distances  along  both  axes  are  in  A. 
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Figure  6.  Intramolecular  interaction  potentials  of  Model  II  for  different  values  (Equation  2) 

as  functions  of  intemuclear  distance  for  (a)  A-B  interactions;  (b)  A-A  interactions; 
and  B-R  interactions.  Magnitude  of  the  values  are  denoted  in  the  legend. 
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The  two-dimensional  crystal  modeled  by  Equation  (2)  also  has  an  equilibrium  configuration  consistent 
with  a  herringbone  arrangemenL  We  determined  the  1  atm,  20  K  orientation  of  the  crystal  using  the  NPT 
Monte  Carlo  calculations  (described  in  Rice  et  al.,  in  press).  The  low-temperature,  ambient-pressure  lattice 
parameters  of  the  crystal  in  the  x  and  y  directions,  respectively,  are  a  =  4.34  A,  b  =  6.27  A.  The  COB 
fractionals  for  the  two  molecules  (1  and  2)  in  the  unit  cell  are  at  0.25,  0.25  and  0.75,  0.75,  respectively. 
The  A-B  equilibrium  bond  length  is  1.349  A,  and  the  molecular  orientation  of  molecules  1  and  2  are  29.0° 
and  -29.0°,  respectively,  relative  to  the  crystal  x-axis. 

NPT  Metropolis  Monte  Carlo  calculations  of  equilibrium  crystal  structure  as  a  function  of  pressure 
at  20  K  were  performed  to  determine  the  sound  speeds  of  the  different  models.  Results  of  the  density  as 
a  function  of  pressure  for  Models  I  and  II  are  shown  in  Table  3.  These  results  for  P  ranging  from  0.0  to 
0.015  eV/A^  were  fitted  to  a  cubic  polynomial,  and  the  slope  of  P-V  curve  was  extracted,  which  is  directly 
proportional  to  the  square  of  the  sound  speed  of  the  crystal  (Pickett  and  Davis  1979;  Pickett  1984).  The 
sound  speeds  for  Models  I  and  n  are  2.0  and  1.2  km/s,  respectively. 

3.  DETAILS  OP  THE  CALCULATIONS 

Details  of  the  molecular  dynamics  simulations  are  given  in  the  accompanying  paper  on  the  comparison 
of  molecular  dynamics  simulations  to  hydrodynamic  predictions  (Rice  et  al.,  in  press).  The  results  that 
are  analyzed  in  this  paper  are  from  the  simulation  in  which  shock  is  initiated  by  a  plate  of  A-A  molecules 
striking  the  simulation  cell  with  an  impact  velocity  of  12  km/s. 

Because  our  interest  is  to  examine  the  mechanistic  details  involved  in  the  reactions  initiating  the 
detonation,  we  will  focus  our  analysis  on  a  region  we  caU  the  reaction  zone.  The  reaction  zone  here  is 
defined  as  the  region  between  the  shock  fiont  and  the  point  at  which  the  number  of  reacted  A-B  molecules 
equals  or  exceeds  the  number  of  unreacted  A-B  molecules.  It  is  difficult  to  determine  whether  a  molecule 
has  reacted  or  not  witiiin  this  zone  in  an  unambiguous  manner,  due  to  the  effect  on  Equation  (2)  by  the 
higher  density  of  this  region.  As  shown  earlier,  there  exists  a  density  of  nuclei  [reflected  by  the  term 

of  Equation  (2)]  that  eliminates  any  attraction  between  the  A-B  pair,  even  though  the  atoms  are  within 
the  3.0  A  cutoff  for  the  intramolecular  term.  However,  we  know  unambiguously  that  some  sort  of 
molecular  interaction  exists  if  the  A-B  intemuclear  distance  is  less  than  3.0  A  according  to  Equation  (2). 
Therefore,  for  simplicity,  we  have  used  the  following  simple  geometric  test  to  determine  reactivity.  We 
have  defined  that  the  original  A-B  pair  has  undergone  reaction  (is  dissociated,  or  is  bound  to  an  atom  that 
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is  not  its  original  molecular  partner)  if  the  distance  between  the  atomic  centers  exceeds  3.0  A.  Figure  8(d) 
of  Rice  et  al.  (in  press)  shows  the  number  of  unreacted  A-B  molecules  and  dissociated  A-B  molecules  at 
7.8  ps  in  the  simulation.  The  point  at  which  the  number  of  dissociated  A-B  molecules  exceeds  the  number 
of  unreacted  A-B  molecules  falls  slightly  behind  the  shock  fronts  in  the  density,  pressure,  and  temperature 
profiles.  As  shown  in  Rice  et  aL  (in  press),  the  width  of  the  reaction  zone  is  approximately  constant 
throughout  die  simulations  of  unsupported  detonations  and  is,  on  average,  14  A  wide. 

We  win  examine  the  translational,  rotational,  and  vibrational  kinetic  energy  (KE)  distributions  of  all 
unreacted  A-B  molecules  within  the  reaction  zone.  Before  the  KE  components  for  each  molecule  are 
calculated,  the  local  mass  flow  velocity  is  removed  as  described  in  Rice  et  al.  (in  press).  Additionally, 
we  will  calculate  the  distributions  of  A-B  bond  lengths  and  orientation  angles  of  the  A-B  molecules 
relative  to  the  crystal  x-axis  in  this  region.  We  will  also  calculate  the  distributions  of  the  intramolecular 
interaction  potential  and  the  term  for  the  unreacted  A-B  molecules.  The  calculation  of  the  translational 

energy  of  the  unreacted  A-B  molecule  is  straightforward: 

p2 

KE  (Translation)  = 

2M 


where  M  is  the  total  mass  of  the  A-B  molecule  and  is  the  center-of-mass  momentum  of  the  molecule. 
As  an  s^proximation,  we  assume  that  the  rotational  energy  is  the  KE  due  to  motion  of  the  atoms  in  the 
molecule  that  are  perpendicular  to  the  molecular  bond 


KE  (Rotation)  =  1  E 
2  i  =  A 


B 


m-v. 

^  ^  (peipendicuiar) 


(4) 


where  Vj  (perpendicular)  denotes  the  velocity  components  perpendicular  to  the  molecular  bond  (center-of-mass 
velocity  removed).  The  vibrational  KE  is  the  KE  due  to  motion  of  the  atoms  in  the  molecule  along  the 
molecular  bond 


KE  (Vibration)  = 


1  ® 

1  E 

2  i  =  A 


m^Vj  , 

1  *  (parallel) 


(5) 
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where  Vj^pajaUei)  denotes  the  velocity  components  parallel  to  the  molecular  bond  (center-of-mass  velocity 
removed). 

Although  these  calculations  of  the  internal  energy  distributions  are  based  on  a  simple  idea,  especially 
in  light  of  the  fact  that  the  molecules  might  not  even  be  bound  in  this  zone  due  to  the  high-density  effects 
on  the  interaction  potential,  it  represents  adequately  energy  flow  into  the  directional  components  of  the 
system. 

4.  RESULTS 

Each  of  the  calculations  that  simulated  unsupported  detonations  discussed  in  Rice  et  al.  (in  press) 
(flyer-plate  impact  velocities  >  4.7  km/s)  provide  the  same  temperature,  density,  and  pressure  profiles  as 
well  as  the  same  microscopic  behavior  in  the  reaction  zone  once  the  detonation  wave  reaches  the  steady- 
state  velocity.  Therefore,  we  are  presenting  the  results  of  only  one  simulation,  the  simulation  with  plate 
impact  of  12  km/s. 

We  have  examined  in  detail  the  various  properties  within  the  reaction  zone  in  order  to  detect  any 
variations  occurring  there.  The  reaction  zone  is  partitioned  into  seven  regions,  labeled  0-6.  Each  region 
has  a  2.17  A  width  in  the  crystal  x  direction  (1/2  the  size  of  the  unit  cell  in  the  x  direction),  and  because 
the  shock  proffle  is  uniform  in  the  y  direction,  the  length  in  the  y-direction  is  the  entire  crystal  length. 
Region  0  is  the  area  directly  behind  the  shock  front  to  2.17  A  behind  the  shock  front;  region  1  begins  at 
the  end  of  region  0,  and  ends  at  4.34  A  behind  the  shock  front;  region  2  ends  6.51  A  behind  the  shock 
front;  and  so  on,  concluding  with  region  6,  which  ranges  from  13.02  to  15.19  A  behind  the  shock  front. 
Distributions  of  properties  of  the  system  were  calculated  at  every  100th  integration  step  throughout  the 
12  km/s  trajectory;  cumulative  distributions  of  these  properties  were  compiled,  normalized,  and  analyzed. 
The  normalization  is  such  that  the  integral  of  the  property  over  the  independent  variable  is  one.  Within 
the  statistical  uncertainty  resulting  from  this  limited  sampling,  none  of  the  properties  changed  in  time  after 
the  first  0.25  ps  following  flyer-plate  impact.  In  order  to  obtain  better  statistical  sampling,  we  therefore 
report  cumulative  normalized  distributions  and  corresponding  time  averages  for  all  times  later  than  0.5  ps 
after  flyer-plate  impact  We  have  calculated  energy  and  orientational  distributions  of  the  unreacted  A-B 
molecule  that  are  in  each  region.  Additionally,  we  have  calculated  the  distributions  of  the  intramolecular 
interactions  and  corresponding  terms  for  the  unreacted  A-B  pairs  in  these  regions.  Averages  obtained 


20 


from  the  normalized  distributions  of  these  properties  for  each  of  the  seven  regions  within  the  reaction  zone 
are  shown  in  Table  4. 

Figures  7-9  show  the  cumulative  translational,  rotational,  and  vibrational  KE  distributions, 
respectively,  of  the  unreacted  A-B  pairs  in  the  seven  zones  throughout  the  12  km/s  simulation.  It  is 
evident  that  the  translational  modes  of  the  A-B  molecules  in  the  region  directly  behind  the  shock  front 
(region  0)  are  substantially  excited  as  the  shock  wave  passes  through  it,  but  the  internal  modes  are  excited 
to  only  one-quarter  of  the  value  of  the  translational  excitation.  The  average  translational  eneigy  of 
molecules  within  region  1  is  less  than  half  that  of  region  0,  while  the  average  energy  of  the  internal  modes 
has  doubled,  indicating  that  substantial  energy  transfer  among  molecular  modes  has  occtirred  by  the  time 
the  molecules  have  moved  from  region  0  to  region  1.  Since  the  average  residence  time  of  the  molecules 
in  each  2.17  A  region  in  the  reaction  zone  is  only  3.3  x  10~^^  s  (assiraiing  local  mass  flow  velocity  of 
6.6  km/s,  the  speed  of  the  detonation  wave  in  this  simulation  [Rice  et  al.,  in  press]),  energy  transfer  among 
molecular  modes  is  rapid.  Beyond  region  1,  the  regional  averages  reach  plateau  values,  consistent  with 
equipartitioning  of  energy  among  translational  and  internal  modes.  It  is  evident  that  substantial  vibrational 
excitation  does  not  occur  due  to  shock  impact.  Therefore,  the  mechanism  of  vibrational  excitation  to 
dissociation  of  A-B  molecules  from  shock  impact  (thermal  decomposition)  is  not  consistent  with  this 
model. 

Figures  10  and  11  show  the  distributions  of  molecular  bonding  energies  for  atom  pairs  and  thelf^ 

terms  wittiin  the  seven  regions.  Rgures  12  and  13  show  the  distributions  of  the  intemuclear  A-B  distances 
and  rotational  orientations  within  the  seven  regions.  The  distribution  of  the  values  in  region  0  is 
spread  across  the  range  of  values,  reflecting  high  density  behind  the  shock  front,  but  has  a  broad 
maximum  at  0.9.  The  average  value  for  region  0  is  0.7.  Rgure  6  shows  that  the  maximum  A-B 

intramolecular  attraction  corresponding  to  this  B^  value  is  -0.3  eV  when  =  1.6  A.  The 
corresponding  distribution  for  intemuclear  distances  of  the  A-B  molecules  within  region  0  (Rgure  12), 
however,  is  a  narrow  distribution  centered  at  1.35  A.  A-B  molecules  with  this  intemuclear  distance  and 
a  B^  value  of  0.7  (or  less)  experience  no  attraction  and  different  degrees  of  repulsion,  depending  on  the 
actual  B^  value.  This  is  reflected  in  the  intramolecular  potential  energy  distribution  for  the  A-B  pairs 
shown  in  Figure  10.  The  distribution  has  a  peak  at  Vjj  =  0  eV  for  region  0,  but  also  has  a  shoulder 
extending  from  V  =  0  to  V  =  -1.  The  shoulder  is  due  to  the  significant  number  of  A-B  molecules  that 
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KE  (Translation)  (eV) 


Figure  7.  Normalized  translational  energy  distributions  of  original  A-B  pairs  for  seven  regions  directly 
behind  the  shock  front  through  the  12  km/s  simulation  (up  to  7.8  ps).  See  text  for  definition 
of  regions. 


KE  (Rotation)  (eV) 


Figure  8.  Normalized  rotational  enerev  distributions  of  original  A-B  pairs  for  seven  regions  directlv 
behind  the  shock  front  through  the  12  km/s  simulation  tun  to  7.8  os).  See  text  for  definition 


of  regions. 
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Figure  9.  Normali/ed  vihrational  KE  distributions  of  original  A-B  pairs  for  seven  regions  directly  behind 
the  shock  front  through  the  12  km/s  simulation  (up  to  7.8  os).  See  text  for  definition  of 
regions. 


Figure  10.  Normalized  distribution  of  intramolecular  terms  of  Equation  (2)  of  original  A-B  pairs  for 
seven  regions  directly  behind  the  shock  front  through  the  12  km/s  simulation  (up  to  7.8  ps). 
See  text  for  definition  of  regions. 
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directlv  behind  the  shock  front  through  the  12  kin/s  simulation  (up  to  7.8  ps).  See  text  for 
definition  of  regions. 


directlv  behind  the  shock  front  through  the  12  km/s  simulation  (up  to  7.8  os’).  See  text  for 
definition  of  regions. 
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7.8  ps).  See  text  for  definition  of  regions. 


have  values  >  0.7  within  region  0.  The  distribution  of  values  within  region  1  is  substantially 
different  fiom  that  of  region  0,  with  its  peak  at  0.1,  but  with  a  broad  shoulder  extending  to  0.9.  The 
average  B^  value  within  this  region  is  0.3,  and  the  corresponding  average  intemuclear  distance  is  1.8  A. 

The  distribution  of  R^bS,  however,  shows  a  peak  at  1.35  A,  with  a  broad  shoulder  extending  to  the 
dissociation  limit  of  3.0  A.  The  intramolecular  attraction  for  a  B^  value  of  0.3  is  slight  (<-0.05  eV  at 
R(eq)  =  2.0  A),  and  the  distribution  of  intemuclear  distances  suggests  that  a  substantial  number  of  A-B 
molecules  experience  no  attraction  and  some  degree  of  repulsion  within  this  zone,  also  reflected  in  the 
distribution  of  the  intramolecular  interactions  in  Figure  10  for  region  1.  Beyond  region  1,  theB^ 
distributions  are  similar,  with  peaks  at  0.05.  Likewise,  the  intemuclear  distance  distributions  are  similar, 
and  are  broad  bands  extending  from  1.2  A  to  the  dissociation  limit  of  3.0  A. 


Figure  13  shows  the  distribution  of  the  orientation  angle  of  the  A-B  molecules  throughout  the  reaction 
zone.  Stracture  in  the  orientational  distributions  exists  in  regions  0-2.  In  the  undisturbed  crystal,  the 
molecular  bond  orientation  is  ±29®  relative  to  the  crystal  x-axis.  The  distribution  of  the  angular 
orientation  in  region  0  is  peaked  at  30°,  indicating,  for  the  most  part,  no  rotational  reorientation  upon 
shock  impact.  This,  coupled  with  the  intemuclear  distance  distribution  for  the  A-B  molecules  (Figure  12) 
and  the  density  of  the  system  for  this  region  (~7  amu/A^)  (Rice  et  al.,  in  press),  indicates  fiiat  the  first 
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geometric  parameters  affected  by  the  shock  wave  are  the  intermolecular  distances  (they  are  reduced).  By 
region  1,  however,  die  angular  orientational  distribution  has  a  broad  peak  at  60°,  showing  that  the  A-B 
molecules  are  lining  up  such  that  the  molecular-bond  axes  are  nearly  perpendicular  to  the  direction  of  the 
shock  wave  propagation.  By  region  3,  the  orientational  distribution  is  flat,  indicating  that  the  original  A-B 
pairs  have  experienced  complete  rotational  disorder.  By  this  region,  atomization  of  the  species  has 
occurred  and  the  atoms  are  energetically  free  to  associate  with  potential  molecular  partners.  We  wish  to 
point  out  that  compression  of  the  herringbone  lattice  in  the  direction  of  shock  propagation,  followed  by 
realignment  of  the  molecular-bond  axes  approximately  perpendicular  to  this  direction  places  the  unreacted 
A-B  molecules  in  an  optimum  arrangement  to  form  the  homonuclear  products.  As  shown  previously,  the 
properties  of  the  intramolecular  interaction  potentials  for  both  heteronuclear  and  htanonuclear  interactions 
are  such  that  at  any  degree  of  density,  homonuclear  attractions  are  greater  than  heteronuclear  and 
homonuclear  repulsions  are  less  than  heteronuclear.  It  is  possible  that  the  compression  followed  by 
realignment  of  the  molecules  initially  behind  the  shock  fiont  might  play  some  role  in  the  rapidity  with 
which  energy  redistribution  among  aU  modes  and  completion  of  reactions  have  occurred. 

These  results  clearly  indicate  that  A-B  molecules  within  the  reaction  zone  are  transitiotting  to  an 
atomic  phase  due  to  the  changes  in  the  intramolecular  interactions  caused  by  the  high  density  behind  the 
shock  front.  From  as  close  at  2.0  A  behind  the  shock  front,  the  A-B  molecules  experience,  on  average, 
repulsive  interactions.  Therefore,  we  have  a  picture  of  A  atoms  and  B  atoms  that  are  in  close  proximity 
to  one  another,  but  are  not  boimd.  For  the  regions  beyond  0  and  1,  the  average  values  and 
corresponding  distributions  show  little  intemuclear  attraction  whatsoever  among  the  atoms,  indicating 
complete  atomization  of  the  products.  As  the  shock  front  proceeds  through  the  crystal,  the  rarefacting 
flow  becomes  less  dense,  increasing  the  values.  Correspondingly,  the  repulsive  intramolecular 

interactions  decrease  while  attractions  increase.  Because  the  well  depths  of  the  homonuclear  species  are 
larger  than  that  of  the  heteronuclear  molecules,  homonuclear  product  formation  is  energetically  favored; 
the  realignment  of  the  molecules  in  region  1  puts  the  system  in  optimum  orientation  for  homonuclear 
product  formation.  All  of  these  factors  result  in  complete  reaction  of  the  shocked  A-B  molecular  crystal 
to  exothermically  form  the  homonuclear  diatomic  products. 

Atomization  of  a  molecular  solid  subjected  to  hydrostatic  high  pressure  is  not  unknown.  Such 
phenomena  have  been  observed  in  high-pressure  studies  of  atomization  of  solid  iodine  (Takemura  et  al. 
1980,  1982;  Fujii  et  al.  1987)  HI  (van  Straaten  and  Silvera  1986),  and  IBr  (Fujii  1985).  Althou^ 
pressure-induced  atomization  was  not  specifically  investigated  in  a  shock-loading  experiment  of  iodine 
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(Alder  and  Christian  1960;  McMahan,  Hoid,  and  Ross  1977),  x-ray  spectra  taken  from  these  experiments 
are  consistent  with  the  x-ray  spectra  of  the  atomized  phase  of  solid  iodine  at  pressures  above  21  GPa 
Takemura  et  al.  1980).  Ab  initio  pseudopotential  calculations  have  shown  pressure-induced  changes  in 
electronic  structure  of  molectilar  crystals,  including  loss  of  covalency  leading  to  molecular  dissociation 
with  increasing  pressure  (McMahan  1986).  Electron  densities  of  these  materials  (Br2, 12.  IBr)  determined 
from  x-ray  diffraction  studies  as  a  function  of  pressure  also  show  delocalization  of  the  covalent  electrons 
(Fujii  1995).  These  rearrangements  of  the  electrons  are  consistent  with  the  present  study  and  provide  the 
mechanism  for  the  transition  from  a  molecular  solid  to  the  atomized  solid  in  several  halogen  materials. 
Experimental  studies  of  energetic  materials  subjected  to  high  pressure  have  shown  that  pressure  can 
drastically  affect  thermal  decomposition  rates  and  reaction  mechanisms  (Piermarini,  Block,  and  Miller 
1987,  1989).  For  example,  the  activation  energy  for  the  thermal  decomposition  of  P-HMX 
(octahydro-l,3,5,7-tetranitro-l,3,5,7,tetrazocine),  one  of  the  most  widely  used  explosives,  approaches  zero 
at  very  high  pressures  (Piermarini,  Block,  and  Miller  1987).  This  was  attributed  to  changes  in  the 
interaction  potential  (i.e.,  electronic  structure)  of  the  system  under  compression.  Similarly,  the  thermal 
rate  of  decomposition  of  nitromethane,  a  prototypical  ejqjlosive,  and  its  chemical  reactivity  were  enhanced 
with  increased  pressure  (Piermarini,  Block,  and  Miller  1989).  These  studies,  both  experimental  and 
theoretical,  on  real  molecular  crystals  indicate  changes  in  the  electronic  structure  of  molecular  bonds  due 
to  pressure  and  suggest  to  us  that  the  pressure-induced  atomization  mechanism  of  detonation  for  our  model 
energetic  crystal  is  not  uiueasonable.  We  stress,  however,  that  the  interaction  potential  for  our  system  was 
not  developed  or  parameterized  to  yield  measured  high-pressure  behavior.  Rather,  the  reaction  mechanism 
(pressure-induced  atomization)  is  a  result  of  the  functional  form  and  choice  of  parameters  for  the  potential. 
It  may  be  that  this  model  fortuitously  describes  correctly  a  significant  reaction  mechanism  occurring  in 
a  detonation.  The  molecules  used  in  the  present  study  are  very  simple,  containing  but  a  single  degree  of 
vibrational  motion.  It  will  be  interesting  to  investigate,  with  these  techniques,  the  effect  of  large  molecules 
(many  vibrational  modes)  on  the  energy-sharing  process  and  their  importance  on  the  reaction  mechamsms. 
It  is  imperative  that  either  new  experimental  techniques  be  developed  to  probe  microscopically  the  state 
of  the  system  behind  a  detonation  wave  or  first-principles  calculations  such  as  those  done  by  McMahan 
(1986)  must  be  carried  out  to  determine  quantitatively  the  changes  in  the  electronic  structure  of  explosives 
when  subjected  to  pressures  consistent  with  those  measured  in  detonations.  Only  then  will  it  be  possible 
to  state  whether  models  such  as  that  presented  in  this  work  are  adequate  to  represent  correctly  the 
chemical  reactions  leading  to  detonation. 
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5.  CONCLUSIONS 


We  have  presented  a  comparative  study  of  models  that  produce  features  of  an  unsupported  detonation. 
The  three  interaction  potentials,  based  on  Abell  (1985)  and  Tersoff’s  (1988)  ideas  and  developed  by 
Brenner  (1992,  1993),  Robertson  et  al.  (1992),  and  White  et  al.  (1992,  1994),  incorporate  many-body 
effects  on  the  intramolecular  terms.  A  many-body  term  attenuates  the  molecular  attractions  with 
increasing  density  while  increasing  the  repulsive  portion  of  the  term.  This  changes  the  equilibrium 
molecular  size  and  binding  energy  of  a  molecule  upon  density  increases.  When  the  density  becomes  large 
enough,  the  atoms  in  molecule  are  pushed  apart  and  feel  no  attraction.  In  other  words,  this  model  predicts 
pressure-induced  atomization. 

Results  of  a  two-dimensional  molecular-dynamics  simulation  of  a  detonation  in  an  accompanying  study 
show  that  for  flyer  plates  moving  with  velocities  greater  than  4.7  km/s,  unsupported  detonations  are 
sustained,  and  have  similar  macroscopic  properties  once  the  detonation  waves  reach  the  steady-state 
velocity  (Rice  et  al.,  in  press).  The  detonation  wave  propagates  through  the  quiescent  crystal  at  6.6  km/s; 
immediately  behind  the  detonation  wave,  there  is,  on  average,  a  14-A-wide  region  in  which  microscopic 
and  macroscopic  properties  of  the  crystal  are  time  independent.  It  is  within  this  region  that  the 
heteronuclear  molecules  are  undergoing  reaction;  we  have  labeled  this  the  reaction  zone. 

We  examined  in  detail  seven  regions  within  the  14-A-reaction  zone  overtime;  cumulative  distributions 
of  molecular  properties  including  translational  and  internal  KEs,  geometric  parameters,  and  intramolecular 
potential  energies  were  calculated.  The  average  values  of  properties  of  the  reaction  zone  converge  to  time- 
independent,  constant  values  approximately  4  A  behind  the  shock  front,  'i^thin  the  4-A  region 
immediately  behind  the  shock  front,  the  shock  wave  first  compresses  the  crystal  and  rotates  the  molecules 
such  that  the  molecular  bonds  are  almost  perpendicular  to  the  shock  wave  propagation.  Translational 
modes  of  the  molecules  are  excited,  but  the  internal  modes  of  the  molecules  are  cold.  Thereafter,  energy 
rapidly  transfers  into  internal  modes  until  an  equipartitiorting  among  these  kinetic  degrees  of  freedom  are 
established.  By  this  point,  the  atoms  experience  an  almost  completely  repulsive  intramolecular  interaction; 
in  other  words,  they  become  atomized.  As  the  shock  wave  passes,  the  density  in  the  rarefacting  flow 
becomes  lower,  concurrently  the  molecular  bonding  attractions  increase,  resulting  in  the  association  of  the 
atoms  to  diatomic  products.  The  attractive  portions  of  the  homonuclear  interaction  potentials  ate 
considerably  larger  than  the  heteronuclear  interaction  potential;  thus,  homonuclear  product  formation  is 
energetically  favored  over  heteronuclear  recombination.  The  homonuclear  association  reactions  release 


29 


substantial  chemical  energy  which  drives  the  detonation  wave.  Vibrational  excitation  of  the  reactant 
heteronuclear  molecules  to  dissociation  (thermal  decomposition)  does  not  contribute  to  the  detonation 
mechanism.  Rather,  the  mechanism  of  detonation  for  this  model  is  pressure-induced  atomization  of  the 
molecules  behind  the  shock  front,  followed  by  association  of  the  atoms  to  form  homonuclear  products. 
This  mechanism  is  attributed  directly  to  the  functional  form  of  the  many-body  term  of  the  interaction 
potential;  subsequent  first-principles  calculations  must  be  done  on  systems  subjected  to  high  pressures  to 
determine  if  this  interaction  potential  correctly  describes  a  significant  mechanism  for  initiation  reactions 
in  detonation. 
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